Chapter 6 Diversity analysis

6.1 Alpha diversity

# Calculate Hill numbers
richness <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 0) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(richness = 1) %>%
  rownames_to_column(var = "sample")

neutral <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(neutral = 1) %>%
  rownames_to_column(var = "sample")

phylogenetic <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1, tree = genome_tree) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(phylogenetic = 1) %>%
  rownames_to_column(var = "sample")

# Aggregate basal GIFT into elements
dist <- genome_gifts %>%
  to.elements(., GIFT_db) %>%
  traits2dist(., method = "gower")

functional <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1, dist = dist) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(functional = 1) %>%
  rownames_to_column(var = "sample") %>%
  mutate(functional = if_else(is.nan(functional), 1, functional))

# Merge all metrics
alpha_div <- richness %>%
  full_join(neutral, by = join_by(sample == sample)) %>%
  full_join(phylogenetic, by = join_by(sample == sample)) %>%
  full_join(functional, by = join_by(sample == sample))

6.1.1 Wild samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="0_Wild") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = Population, group=Population, color=Population, fill=Population)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#008080', "#d57d2c")) +
      scale_fill_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#00808050', "#d57d2c50")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.58) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

6.1.2 Acclimation samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="1_Acclimation") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = Population, group=Population, color=Population, fill=Population)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#008080', "#d57d2c")) +
      scale_fill_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#00808050', "#d57d2c50")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.58) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

6.1.3 Antibiotics samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="2_Antibiotics") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = Population, group=Population, color=Population, fill=Population)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#008080', "#d57d2c")) +
      scale_fill_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#00808050', "#d57d2c50")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.58) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

6.1.4 Transplant_1 samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="3_Transplant1") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA","#d57d2c", "#76b183")) +
      scale_fill_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA50","#d57d2c50","#76b18350")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.7) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

6.1.5 Transplant_2 samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="4_Transplant2") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA","#d57d2c", "#76b183")) +
      scale_fill_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA50","#d57d2c50","#76b18350")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.7) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

6.1.6 Post-Transplant_1 samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="5_Post-FMT1") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA","#d57d2c", "#76b183")) +
      scale_fill_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA50","#d57d2c50","#76b18350")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.7) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

6.1.7 Post-Transplant_2 samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="6_Post-FMT2") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA","#d57d2c", "#76b183")) +
      scale_fill_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA50","#d57d2c50","#76b18350")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.7) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

6.2 Beta diversity

beta_q0n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 0)

beta_q1n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1)

beta_q1p <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1, tree = genome_tree)

beta_q1f <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1, dist = dist)

6.3 Permanovas

6.3.0.1 Load required data

meta <- column_to_rownames(sample_metadata, "Tube_code")

6.3.1 1. Are the wild populations similar?

6.3.1.1 Wild: P.muralis vs P.liolepis

wild <- meta %>%
  filter(time_point == "0_Wild")

# Create a temporary modified version of genome_counts_filt
temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

wild.counts <- temp_genome_counts[, which(colnames(temp_genome_counts) %in% rownames(wild))]
identical(sort(colnames(wild.counts)), sort(as.character(rownames(wild))))

wild_nmds <- sample_metadata %>%
  filter(time_point == "0_Wild")

6.3.1.2 Number of samples used

[1] 27

beta_div_richness_wild<-hillpair(data=wild.counts, q=0)
beta_div_neutral_wild<-hillpair(data=wild.counts, q=1)
beta_div_phylo_wild<-hillpair(data=wild.counts, q=1, tree=genome_tree)
beta_div_func_wild<-hillpair(data=wild.counts, q=1, dist=dist)

6.3.1.3 Richness

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.000012 0.000012 0.0012    999   0.98
Residuals 25 0.257281 0.010291                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
                  Podarcis_liolepis Podarcis_muralis
Podarcis_liolepis                               0.98
Podarcis_muralis            0.97302                 

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	1.542719	0.2095041	6.625717	0.001
Residual	25	5.820951	0.7904959	NA	NA
Total	26	7.363669	1.0000000	NA	NA

6.3.1.4 Neutral

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.000048 0.0000476 0.0044    999   0.94
Residuals 25 0.270114 0.0108046                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
                  Podarcis_liolepis Podarcis_muralis
Podarcis_liolepis                               0.94
Podarcis_muralis            0.94763                 

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	1.918266	0.2608511	8.822682	0.001
Residual	25	5.435610	0.7391489	NA	NA
Total	26	7.353876	1.0000000	NA	NA

6.3.1.5 Phylogenetic

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.03585 0.035847 2.4912    999  0.127
Residuals 25 0.35973 0.014389                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
                  Podarcis_liolepis Podarcis_muralis
Podarcis_liolepis                              0.127
Podarcis_muralis            0.12705                 

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	0.3218613	0.2162815	6.899207	0.001
Residual	25	1.1662981	0.7837185	NA	NA
Total	26	1.4881594	1.0000000	NA	NA

6.3.1.6 Functional

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.018367 0.018367 1.5597    999  0.245
Residuals 25 0.294402 0.011776                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
                  Podarcis_liolepis Podarcis_muralis
Podarcis_liolepis                              0.246
Podarcis_muralis            0.22328                 

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	0.0858578	0.172879	5.225323	0.045
Residual	25	0.4107775	0.827121	NA	NA
Total	26	0.4966352	1.000000	NA	NA

beta_q0n_nmds_wild <- beta_div_richness_wild$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(wild_nmds, by = join_by(sample == Tube_code))

beta_q1n_nmds_wild <- beta_div_neutral_wild$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(wild_nmds, by = join_by(sample == Tube_code))

beta_q1p_nmds_wild <- beta_div_phylo_wild$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(wild_nmds, by = join_by(sample == Tube_code))

beta_q1f_nmds_wild <- beta_div_func_wild$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(wild_nmds, by = join_by(sample == Tube_code))

6.3.2 2. Effect of acclimation

accli <- meta %>%
  filter(time_point == "1_Acclimation")

# Create a temporary modified version of genome_counts_filt
temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

accli.counts <- temp_genome_counts[, which(colnames(temp_genome_counts) %in% rownames(accli))]
identical(sort(colnames(accli.counts)), sort(as.character(rownames(accli))))

accli_nmds <- sample_metadata %>%
  filter(time_point == "1_Acclimation")

6.3.2.1 Number of samples used

[1] 27

beta_div_richness_accli<-hillpair(data=accli.counts, q=0)
beta_div_neutral_accli<-hillpair(data=accli.counts, q=1)
beta_div_phylo_accli<-hillpair(data=accli.counts, q=1, tree=genome_tree)
beta_div_func_accli<-hillpair(data=accli.counts, q=1, dist=dist)

6.3.2.2 Richness

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)   
Groups     1 0.11796 0.117959 12.963    999  0.002 **
Residuals 25 0.22748 0.009099                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
          Cold_wet Hot_dry
Cold_wet             0.002
Hot_dry  0.0013711        

	Df	SumOfSqs	R2	F	Pr(>F)
Population	1	1.639807	0.179834	5.481634	0.001
Residual	25	7.478640	0.820166	NA	NA
Total	26	9.118447	1.000000	NA	NA

6.3.2.3 Neutral

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.07844 0.078443 5.2384    999  0.025 *
Residuals 25 0.37437 0.014975                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.029
Hot_dry  0.030815        

	Df	SumOfSqs	R2	F	Pr(>F)
Population	1	1.947003	0.2306127	7.493387	0.001
Residual	25	6.495736	0.7693873	NA	NA
Total	26	8.442739	1.0000000	NA	NA

6.3.2.4 Phylogenetic

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.06739 0.067395 2.9532    999  0.102
Residuals 25 0.57052 0.022821                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.101
Hot_dry  0.098068        

	Df	SumOfSqs	R2	F	Pr(>F)
Population	1	0.2441653	0.1224638	3.488854	0.03
Residual	25	1.7496100	0.8775362	NA	NA
Total	26	1.9937754	1.0000000	NA	NA

6.3.2.5 Functional

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.02496 0.024955 0.6729    999   0.42
Residuals 25 0.92714 0.037085                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.413
Hot_dry   0.41979        

	Df	SumOfSqs	R2	F	Pr(>F)
Population	1	0.0279454	0.0248037	0.6358634	0.425
Residual	25	1.0987171	0.9751963	NA	NA
Total	26	1.1266624	1.0000000	NA	NA

beta_q0n_nmds_accli <- beta_div_richness_accli$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(accli_nmds, by = join_by(sample == Tube_code))

beta_q1n_nmds_accli <- beta_div_neutral_accli$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(accli_nmds, by = join_by(sample == Tube_code))

beta_q1p_nmds_accli <- beta_div_phylo_accli$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(accli_nmds, by = join_by(sample == Tube_code))

beta_q1f_nmds_accli <- beta_div_func_accli$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(accli_nmds, by = join_by(sample == Tube_code))

6.3.3 3. Comparison between Wild and Acclimation

accli1 <- meta  %>%
  filter(time_point == "0_Wild" | time_point == "1_Acclimation")

temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

accli1.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(accli1))]
identical(sort(colnames(accli1.counts)),sort(as.character(rownames(accli1))))

accli1_nmds <- sample_metadata %>%
  filter(time_point == "0_Wild" | time_point == "1_Acclimation")

6.3.3.1 Number of samples used

[1] 54

beta_div_richness_accli1<-hillpair(data=accli1.counts, q=0)
beta_div_neutral_accli1<-hillpair(data=accli1.counts, q=1)
beta_div_phylo_accli1<-hillpair(data=accli1.counts, q=1, tree=genome_tree)
beta_div_func_accli1<-hillpair(data=accli1.counts, q=1, dist=dist)

6.3.3.1.1 Richness

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.05014 0.050145 6.2252    999   0.02 *
Residuals 52 0.41886 0.008055                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
                0_Wild 1_Acclimation
0_Wild                         0.015
1_Acclimation 0.015808              

	Df	SumOfSqs	R2	F	Pr(>F)
time_point	1	0.6172653	0.0360987	2.320618	0.001
species	1	2.8279677	0.1653842	10.631785	0.003
time_point:species	1	0.3545578	0.0207351	1.332965	0.063
Residual	50	13.2995905	0.7777820	NA	NA
Total	53	17.0993812	1.0000000	NA	NA

6.3.3.1.2 Neutral

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.0199 0.0199035 2.1213    999  0.143
Residuals 52 0.4879 0.0093827                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
               0_Wild 1_Acclimation
0_Wild                        0.145
1_Acclimation 0.15128              

	Df	SumOfSqs	R2	F	Pr(>F)
time_point	1	0.9050519	0.0541893	3.792749	0.001
species	1	3.3236300	0.1989999	13.928143	0.001
time_point:species	1	0.5416391	0.0324302	2.269815	0.006
Residual	50	11.9313461	0.7143805	NA	NA
Total	53	16.7016671	1.0000000	NA	NA

6.3.3.1.3 Phylogenetic

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.01334 0.013340 0.6524    999  0.446
Residuals 52 1.06332 0.020449                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
               0_Wild 1_Acclimation
0_Wild                        0.451
1_Acclimation 0.42294              

	Df	SumOfSqs	R2	F	Pr(>F)
time_point	1	0.2890434	0.0766494	4.956318	0.003
species	1	0.3508889	0.0930498	6.016803	0.001
time_point:species	1	0.2151377	0.0570509	3.689034	0.005
Residual	50	2.9159082	0.7732498	NA	NA
Total	53	3.7709782	1.0000000	NA	NA

6.3.3.1.4 Functional

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.0123 0.012300 0.4817    999  0.455
Residuals 52 1.3277 0.025533                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
               0_Wild 1_Acclimation
0_Wild                         0.46
1_Acclimation 0.49073              

	Df	SumOfSqs	R2	F	Pr(>F)
time_point	1	0.0448774	0.0269021	1.4865056	0.236
species	1	0.0973005	0.0583275	3.2229509	0.314
time_point:species	1	0.0165026	0.0098926	0.5466273	0.384
Residual	50	1.5094945	0.9048777	NA	NA
Total	53	1.6681751	1.0000000	NA	NA

beta_richness_nmds_accli1 <- beta_div_richness_accli1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(accli1_nmds, by = c("sample" = "Tube_code"))

beta_neutral_nmds_accli1 <- beta_div_neutral_accli1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(accli1_nmds, by = c("sample" = "Tube_code"))

beta_phylo_nmds_accli1 <- beta_div_phylo_accli1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(accli1_nmds, by = join_by(sample == Tube_code))

beta_func_nmds_accli1 <- beta_div_func_accli1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(accli1_nmds, by = join_by(sample == Tube_code))

6.3.4 4. Do the antibiotics work?

6.3.4.1 Antibiotics

treat1 <- meta  %>%
  filter(time_point == "2_Antibiotics")

temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

treat1.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(treat1))]
identical(sort(colnames(treat1.counts)),sort(as.character(rownames(treat1))))

treat1_nmds <- sample_metadata %>%
  filter(time_point == "2_Antibiotics")

6.3.4.2 Number of samples used

[1] 23

beta_div_richness_treat1<-hillpair(data=treat1.counts, q=0)
beta_div_neutral_treat1<-hillpair(data=treat1.counts, q=1)
beta_div_phylo_treat1<-hillpair(data=treat1.counts, q=1, tree=genome_tree)
beta_div_func_treat1<-hillpair(data=treat1.counts, q=1, dist=dist)

6.3.4.2.1 Richness

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)   
Groups     1 0.015319 0.0153186 6.8764    999   0.01 **
Residuals 21 0.046782 0.0022277                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.011
Hot_dry  0.015919        

	Df	SumOfSqs	R2	F	Pr(>F)
Population	1	1.356644	0.1527052	3.784762	0.001
Residual	21	7.527429	0.8472948	NA	NA
Total	22	8.884073	1.0000000	NA	NA

6.3.4.2.2 Neutral

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.030536 0.0305358 3.8593    999   0.07 .
Residuals 21 0.166158 0.0079123                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.081
Hot_dry  0.062842        

	Df	SumOfSqs	R2	F	Pr(>F)
Population	1	1.785669	0.2085055	5.532084	0.001
Residual	21	6.778468	0.7914945	NA	NA
Total	22	8.564137	1.0000000	NA	NA

6.3.4.2.3 Phylogenetic

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.012041 0.012041 0.9898    999  0.335
Residuals 21 0.255459 0.012165                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.344
Hot_dry   0.33111        

	Df	SumOfSqs	R2	F	Pr(>F)
Population	1	0.8963254	0.1888758	4.889993	0.001
Residual	21	3.8492558	0.8111242	NA	NA
Total	22	4.7455811	1.0000000	NA	NA

6.3.4.2.4 Functional

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.01802 0.018021 0.4386    999  0.524
Residuals 21 0.86280 0.041086                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.537
Hot_dry   0.51499        

	Df	SumOfSqs	R2	F	Pr(>F)
Population	1	0.0184663	0.0098404	0.2087022	0.689
Residual	21	1.8581156	0.9901596	NA	NA
Total	22	1.8765819	1.0000000	NA	NA

beta_richness_nmds_treat1 <- beta_div_richness_treat1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat1_nmds, by = c("sample" = "Tube_code"))

beta_neutral_nmds_treat1 <- beta_div_neutral_treat1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat1_nmds, by = c("sample" = "Tube_code"))

beta_phylo_nmds_treat1 <- beta_div_phylo_treat1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat1_nmds, by = join_by(sample == Tube_code))

beta_func_nmds_treat1 <- beta_div_func_treat1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat1_nmds, by = join_by(sample == Tube_code))

6.3.4.3 Acclimation vs antibiotics

treat <- meta  %>%
  filter(time_point == "1_Acclimation" | time_point == "2_Antibiotics")

temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

treat.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(treat))]
identical(sort(colnames(treat.counts)),sort(as.character(rownames(treat))))

treat_nmds <- sample_metadata %>%
  filter(time_point == "1_Acclimation" | time_point == "2_Antibiotics")

6.3.4.4 Number of samples used

[1] 50

beta_div_richness_treat<-hillpair(data=treat.counts, q=0)
beta_div_neutral_treat<-hillpair(data=treat.counts, q=1)
beta_div_phylo_treat<-hillpair(data=treat.counts, q=1, tree=genome_tree)
beta_div_func_treat<-hillpair(data=treat.counts, q=1, dist=dist)

6.3.4.4.1 Richness

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq     F N.Perm Pr(>F)  
Groups     1 0.025318 0.0253178 6.021    999  0.019 *
Residuals 48 0.201837 0.0042049                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 2_Antibiotics
1_Acclimation                       0.017
2_Antibiotics      0.017817              

	Df	SumOfSqs	R2	F	Pr(>F)
time_point	1	1.8885838	0.0949462	5.789315	0.001
Population	1	2.1171094	0.1064350	6.489843	0.001
time_point:Population	1	0.8793415	0.0442078	2.695557	0.003
Residual	46	15.0060684	0.7544111	NA	NA
Total	49	19.8911031	1.0000000	NA	NA

6.3.4.4.2 Neutral

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.039587 0.039587 6.8387    999  0.015 *
Residuals 48 0.277854 0.005789                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 2_Antibiotics
1_Acclimation                       0.015
2_Antibiotics      0.011886              

	Df	SumOfSqs	R2	F	Pr(>F)
time_point	1	2.0241808	0.1063620	7.014531	0.001
Population	1	2.8531033	0.1499183	9.887052	0.001
time_point:Population	1	0.8795688	0.0462175	3.048030	0.004
Residual	46	13.2742044	0.6975022	NA	NA
Total	49	19.0310573	1.0000000	NA	NA

6.3.4.4.3 Phylogenetic

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq Mean Sq      F N.Perm Pr(>F)    
Groups     1 0.58372 0.58372 35.413    999  0.001 ***
Residuals 48 0.79119 0.01648                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 2_Antibiotics
1_Acclimation                       0.001
2_Antibiotics    2.9795e-07              

	Df	SumOfSqs	R2	F	Pr(>F)
time_point	1	1.8065206	0.2113909	14.842282	0.001
Population	1	0.7903334	0.0924813	6.493340	0.001
time_point:Population	1	0.3501572	0.0409738	2.876874	0.028
Residual	46	5.5988658	0.6551540	NA	NA
Total	49	8.5458771	1.0000000	NA	NA

6.3.4.4.4 Functional

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.18591 0.185914 5.0679    999   0.04 *
Residuals 48 1.76088 0.036685                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 2_Antibiotics
1_Acclimation                       0.038
2_Antibiotics      0.028989              

	Df	SumOfSqs	R2	F	Pr(>F)
time_point	1	1.8020952	0.3750193	28.0355332	0.001
Population	1	0.0031247	0.0006503	0.0486115	0.001
time_point:Population	1	0.0432870	0.0090081	0.6734238	0.498
Residual	46	2.9568327	0.6153223	NA	NA
Total	49	4.8053396	1.0000000	NA	NA

beta_richness_nmds_treat <- beta_div_richness_treat$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat_nmds, by = c("sample" = "Tube_code"))

beta_neutral_nmds_treat <- beta_div_neutral_treat$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat_nmds, by = c("sample" = "Tube_code"))

beta_phylo_nmds_treat <- beta_div_phylo_treat$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat_nmds, by = join_by(sample == Tube_code))

beta_func_nmds_treat <- beta_div_func_treat$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat_nmds, by = join_by(sample == Tube_code))

6.3.5 5. Does the FMT work?

6.3.5.1 Comparison between FMT2 vs Post-FMT2

#Create newID to identify duplicated samples
transplants_metadata<-sample_metadata%>%
  mutate(Tube_code=str_remove_all(Tube_code, "_a"))
transplants_metadata$newID <- paste(transplants_metadata$Tube_code, "_", transplants_metadata$individual)

transplant3<-transplants_metadata%>%
  filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2")%>%
  column_to_rownames("newID")

transplant3_nmds <- transplants_metadata %>%
  filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2")

full_counts<-temp_genome_counts %>%
    t()%>%
    as.data.frame()%>%
    rownames_to_column("Tube_code")%>%
    full_join(transplants_metadata,by = join_by(Tube_code == Tube_code))

transplant3_counts<-full_counts %>%
  filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2") %>%
  subset(select=-c(315:324)) %>%
  column_to_rownames("newID")%>%
  subset(select=-c(1))%>%
  t() %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)

identical(sort(colnames(transplant3_counts)),sort(as.character(rownames(transplant3))))

6.3.5.2 Number of samples used

[1] 49

beta_div_richness_transplant3<-hillpair(data=transplant3_counts, q=0)
beta_div_neutral_transplant3<-hillpair(data=transplant3_counts, q=1)
beta_div_phylo_transplant3<-hillpair(data=transplant3_counts, q=1, tree=genome_tree)
beta_div_func_transplant3<-hillpair(data=transplant3_counts, q=1, dist=dist)

#Arrange of metadata dataframe
transplant3_arrange<-transplant3[labels(beta_div_neutral_transplant3$S),]

6.3.5.2.1 Richness

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	1.180473	0.0855095	5.155229	0.035
time_point	1	0.860906	0.0623612	3.759652	0.001
type	1	1.459433	0.1057165	6.373471	0.005
Residual	45	10.304350	0.7464128	NA	NA
Total	48	13.805162	1.0000000	NA	NA

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control vs Treatment	1	1.4169018	5.739828	0.15622903	0.001	0.003	*
Control vs Hot_control	1	2.0940966	8.509112	0.21005427	0.001	0.003	*
Treatment vs Hot_control	1	0.3004618	1.265034	0.04179854	0.152	0.456

6.3.5.2.2 Neutral

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	1.2800927	0.0939787	6.068484	0.006
time_point	1	0.9350566	0.0686477	4.432785	0.001
type	1	1.9135997	0.1404879	9.071725	0.001
Residual	45	9.4923500	0.6968858	NA	NA
Total	48	13.6210990	1.0000000	NA	NA

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control vs Treatment	1	1.8758788	8.282671	0.21084796	0.001	0.003	*
Control vs Hot_control	1	2.4396317	10.635546	0.24945256	0.001	0.003	*
Treatment vs Hot_control	1	0.3158428	1.394345	0.04587515	0.132	0.396

6.3.5.2.3 Phylogenetic

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	0.1400466	0.0952654	5.873615	0.050
time_point	1	0.1138047	0.0774145	4.773017	0.002
type	1	0.1432667	0.0974558	6.008665	0.004
Residual	45	1.0729504	0.7298643	NA	NA
Total	48	1.4700683	1.0000000	NA	NA

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control vs Treatment	1	0.14387705	5.735321	0.15612552	0.001	0.003	*
Control vs Hot_control	1	0.22715701	9.044894	0.22036587	0.001	0.003	*
Treatment vs Hot_control	1	0.04648319	1.704277	0.05550617	0.152	0.456

6.3.5.2.4 Functional

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	0.0092808	0.0077189	0.3741811	0.396
time_point	1	-0.0061674	-0.0051295	-0.2486581	0.896
type	1	0.0831052	0.0691191	3.3506286	0.226
Residual	45	1.1161295	0.9282915	NA	NA
Total	48	1.2023481	1.0000000	NA	NA

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control vs Treatment	1	0.078539743	4.59293783	0.129040706	0.089	0.267
Control vs Hot_control	1	0.052468954	2.13675422	0.062593948	0.191	0.573
Treatment vs Hot_control	1	-0.002340352	-0.07432315	-0.002569452	0.861	1.000

beta_richness_nmds_transplant3 <- beta_div_richness_transplant3$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(transplant3_nmds, by = join_by(sample == newID))

beta_neutral_nmds_transplant3 <- beta_div_neutral_transplant3$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(transplant3_nmds, by = join_by(sample == newID))

beta_phylo_nmds_transplant3 <- beta_div_phylo_transplant3$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(transplant3_nmds, by = join_by(sample == newID))

beta_func_nmds_transplant3 <- beta_div_func_transplant3$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(transplant3_nmds, by = join_by(sample == newID))

p0<-beta_richness_nmds_transplant3 %>%
            group_by(individual) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p1<-beta_neutral_nmds_transplant3 %>%
            group_by(individual) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")
  
p2<-beta_phylo_nmds_transplant3 %>%
            group_by(individual) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p3<-beta_func_nmds_transplant3 %>%
            group_by(individual) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
                theme_classic()+
                theme(legend.position="none")

6.3.5.3 Comparison between the different experimental time points (Acclimation vs Transplant samples)

The estimated time for calculating the 5151 pairwise combinations is 27 seconds.

ggarrange( p1, p2, p3, ncol=3, nrow=1, common.legend = TRUE, legend="right")

6.3.5.4 Comparison of acclimation samples to transplant samples

transplant7<-transplants_metadata%>%
  filter(time_point == "4_Transplant2" | time_point == "1_Acclimation"| time_point == "3_Transplant1")%>%
  column_to_rownames("newID")

transplant7_nmds <- transplants_metadata %>%
  filter(time_point == "4_Transplant2" | time_point == "1_Acclimation"| time_point == "3_Transplant1")

transplant7_counts<-full_counts %>%
  filter(time_point == "4_Transplant2" | time_point == "1_Acclimation"| time_point == "3_Transplant1") %>%
  subset(select=-c(315:324)) %>%
  column_to_rownames("newID")%>%
  subset(select=-c(1))%>%
  t() %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)

transplant7_counts <- transplant7_counts[, !names(transplant7_counts) %in% c("AD45 _ LI1_2nd_2", "AD48 _ LI1_2nd_6")]

identical(sort(colnames(transplant7_counts)),sort(as.character(rownames(transplant7))))

[1] TRUE

6.3.5.5 Number of samples used

[1] 73

beta_div_richness_transplant7<-hillpair(data=transplant7_counts, q=0)
beta_div_neutral_transplant7<-hillpair(data=transplant7_counts, q=1)
beta_div_phylo_transplant7<-hillpair(data=transplant7_counts, q=1, tree=genome_tree)
beta_div_func_transplant7<-hillpair(data=transplant7_counts, q=1, dist=dist)

#Arrange of metadata dataframe
transplant7_arrange<-transplant7[labels(beta_div_neutral_transplant7$S),]
transplant7_arrange <- transplant7_arrange %>%
  mutate(time_point = recode(time_point,
                             "3_Transplant1" = "Transplant",
                             "4_Transplant2" = "Transplant"))

transplant7_arrange$type_time <- interaction(transplant7_arrange$type, transplant7_arrange$time_point)

6.3.5.5.1 Richness

	Df	SumOfSqs	R2	F	Pr(>F)
Population	1	2.3108184	0.1096208	9.6707812	0.089
time_point	2	1.1082036	0.0525710	2.3189174	0.001
type	1	1.4676332	0.0696217	6.1420489	0.036
Population:time_point	2	0.4228641	0.0200599	0.8848438	0.482
Residual	66	15.7705995	0.7481267	NA	NA
Total	72	21.0801189	1.0000000	NA	NA

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	0.36208146	1.0521088	0.06169963	0.325	1.000
Control.1_Acclimation vs Hot_control.1_Acclimation	1	1.28008774	4.6054436	0.22350616	0.001	0.015	.
Control.1_Acclimation vs Control.Transplant	1	0.55038651	2.2107376	0.08124505	0.003	0.045	.
Control.1_Acclimation vs Treatment.Transplant	1	1.62289430	6.7106689	0.25123553	0.001	0.015	.
Control.1_Acclimation vs Hot_control.Transplant	1	1.73215888	7.4315069	0.25250175	0.001	0.015	.
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	1.36066298	5.0871520	0.24124415	0.001	0.015	.
Treatment.1_Acclimation vs Control.Transplant	1	0.52860586	2.1820402	0.08027507	0.003	0.045	.
Treatment.1_Acclimation vs Treatment.Transplant	1	1.76810026	7.5736721	0.27467042	0.001	0.015	.
Treatment.1_Acclimation vs Hot_control.Transplant	1	1.87790626	8.3291875	0.27462613	0.001	0.015	.
Hot_control.1_Acclimation vs Control.Transplant	1	1.75314247	8.7706781	0.25971282	0.001	0.015	.
Hot_control.1_Acclimation vs Treatment.Transplant	1	0.27700454	1.5346880	0.07126586	0.094	1.000
Hot_control.1_Acclimation vs Hot_control.Transplant	1	0.26448976	1.4916174	0.06349573	0.101	1.000
Control.Transplant vs Treatment.Transplant	1	2.30884687	12.4299510	0.30002331	0.001	0.015	.
Control.Transplant vs Hot_control.Transplant	1	2.50396161	13.6713271	0.30604256	0.001	0.015	.
Treatment.Transplant vs Hot_control.Transplant	1	0.01688622	0.1023282	0.00392027	1.000	1.000

6.3.5.5.2 Neutral

	Df	SumOfSqs	R2	F	Pr(>F)
Population	1	2.641095	0.1266245	12.453113	0.041
time_point	2	1.394944	0.0668791	3.288673	0.001
type	1	1.523725	0.0730534	7.184566	0.031
time_point:type	4	1.724614	0.0826848	2.032946	0.015
Residual	64	13.573319	0.6507583	NA	NA
Total	72	20.857698	1.0000000	NA	NA

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	0.23160196	0.7712905	0.045988741	0.716	1.000
Control.1_Acclimation vs Hot_control.1_Acclimation	1	1.40153474	5.7562378	0.264578733	0.001	0.015	.
Control.1_Acclimation vs Control.Transplant	1	0.56111203	2.5583085	0.092832565	0.001	0.015	.
Control.1_Acclimation vs Treatment.Transplant	1	1.88709838	8.3257794	0.293929402	0.001	0.015	.
Control.1_Acclimation vs Hot_control.Transplant	1	2.02585000	9.2317432	0.295588471	0.001	0.015	.
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	1.63477039	6.8326887	0.299250291	0.001	0.015	.
Treatment.1_Acclimation vs Control.Transplant	1	0.61335323	2.8313912	0.101733730	0.001	0.015	.
Treatment.1_Acclimation vs Treatment.Transplant	1	2.10939140	9.4473664	0.320822116	0.001	0.015	.
Treatment.1_Acclimation vs Hot_control.Transplant	1	2.24827218	10.3907678	0.320794118	0.001	0.015	.
Hot_control.1_Acclimation vs Control.Transplant	1	1.87351542	10.3925002	0.293635661	0.001	0.015	.
Hot_control.1_Acclimation vs Treatment.Transplant	1	0.34276062	1.9273510	0.087897118	0.055	0.825
Hot_control.1_Acclimation vs Hot_control.Transplant	1	0.31638309	1.8072337	0.075911118	0.076	1.000
Control.Transplant vs Treatment.Transplant	1	2.48701901	14.0199769	0.325894571	0.001	0.015	.
Control.Transplant vs Hot_control.Transplant	1	2.75304261	15.6912860	0.336064549	0.001	0.015	.
Treatment.Transplant vs Hot_control.Transplant	1	0.01764676	0.1022118	0.003915827	0.996	1.000

6.3.5.5.3 Phylogenetic

	Df	SumOfSqs	R2	F	Pr(>F)
Population	1	0.3001462	0.0764583	6.402286	0.350
time_point	2	0.2984483	0.0760258	3.183035	0.001
type	1	0.1391084	0.0354361	2.967261	0.626
Residual	68	3.1879143	0.8120798	NA	NA
Total	72	3.9256172	1.0000000	NA	NA

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	0.04186923	0.43916424	0.026714511	0.736	1.000
Control.1_Acclimation vs Hot_control.1_Acclimation	1	0.15609416	2.55468892	0.137684276	0.036	0.540
Control.1_Acclimation vs Control.Transplant	1	0.03888650	0.83961027	0.032493148	0.459	1.000
Control.1_Acclimation vs Treatment.Transplant	1	0.28946588	4.58406811	0.186464994	0.001	0.015	.
Control.1_Acclimation vs Hot_control.Transplant	1	0.31864880	5.37781508	0.196429666	0.001	0.015	.
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	0.23108846	4.05218385	0.202081922	0.003	0.045	.
Treatment.1_Acclimation vs Control.Transplant	1	0.11794420	2.69844074	0.097422117	0.036	0.540
Treatment.1_Acclimation vs Treatment.Transplant	1	0.37640156	6.28511923	0.239113210	0.001	0.015	.
Treatment.1_Acclimation vs Hot_control.Transplant	1	0.40433696	7.18306079	0.246138020	0.001	0.015	.
Hot_control.1_Acclimation vs Control.Transplant	1	0.11597038	5.32063275	0.175478948	0.001	0.015	.
Hot_control.1_Acclimation vs Treatment.Transplant	1	0.03673004	1.13023077	0.053488804	0.330	1.000
Hot_control.1_Acclimation vs Hot_control.Transplant	1	0.04097680	1.30539166	0.056012432	0.240	1.000
Control.Transplant vs Treatment.Transplant	1	0.21736741	7.59281199	0.207494630	0.001	0.015	.
Control.Transplant vs Hot_control.Transplant	1	0.25837791	9.19762187	0.228810100	0.001	0.015	.
Treatment.Transplant vs Hot_control.Transplant	1	0.00180330	0.04804393	0.001844435	0.969	1.000

6.3.5.5.4 Functional

	Df	SumOfSqs	R2	F	Pr(>F)
Population	1	0.0968624	0.0459940	4.124655	0.398
time_point	2	0.1660358	0.0788403	3.535121	0.032
type	1	0.2461830	0.1168973	10.483121	0.182
Residual	68	1.5968952	0.7582683	NA	NA
Total	72	2.1059764	1.0000000	NA	NA

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	0.0904501448	1.65575459	0.0937798825	0.228	1.000
Control.1_Acclimation vs Hot_control.1_Acclimation	1	0.0857692246	1.63605364	0.0927675587	0.225	1.000
Control.1_Acclimation vs Control.Transplant	1	0.0705970944	2.02878350	0.0750601115	0.174	1.000
Control.1_Acclimation vs Treatment.Transplant	1	0.3221196060	7.21718333	0.2651701040	0.006	0.090
Control.1_Acclimation vs Hot_control.Transplant	1	0.3443779759	8.48534465	0.2783417653	0.011	0.165
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	0.0013738865	0.07238159	0.0045034760	0.620	1.000
Treatment.1_Acclimation vs Control.Transplant	1	-0.0041789330	-0.31199952	-0.0126376990	0.779	1.000
Treatment.1_Acclimation vs Treatment.Transplant	1	0.0789743467	4.41748338	0.1809147697	0.085	1.000
Treatment.1_Acclimation vs Hot_control.Transplant	1	0.0842190618	5.17868542	0.1905421598	0.057	0.855
Hot_control.1_Acclimation vs Control.Transplant	1	-0.0042784688	-0.35701715	-0.0144875787	0.818	1.000
Hot_control.1_Acclimation vs Treatment.Transplant	1	0.0825042562	5.11970349	0.2038122581	0.068	1.000
Hot_control.1_Acclimation vs Hot_control.Transplant	1	0.0889170506	6.06518379	0.2161106028	0.040	0.600
Control.Transplant vs Treatment.Transplant	1	0.1900996563	15.59722495	0.3497353247	0.006	0.090
Control.Transplant vs Hot_control.Transplant	1	0.2049565481	17.96453178	0.3668886666	0.004	0.060
Treatment.Transplant vs Hot_control.Transplant	1	-0.0002132009	-0.01472772	-0.0005667719	0.693	1.000

beta_richness_nmds_transplant7 <- beta_div_richness_transplant7$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(transplant7_nmds, by = join_by(sample == newID))

beta_neutral_nmds_transplant7 <- beta_div_neutral_transplant7$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(transplant7_nmds, by = join_by(sample == newID))

beta_phylo_nmds_transplant7 <- beta_div_phylo_transplant7$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(transplant7_nmds, by = join_by(sample == newID))

beta_func_nmds_transplant7 <- beta_div_func_transplant7$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(transplant7_nmds, by = join_by(sample == newID))

p0<-beta_richness_nmds_transplant7 %>%
  group_by(type) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
  scale_color_manual(name="Type",
                     breaks=c("Control", "Hot_control", "Treatment"),
                     labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                     values=c("#4477AA","#d57d2c","#76b183")) +
  scale_shape_manual(name="time_point",
                     breaks=c("1_Acclimation", "3_Transplant1", "4_Transplant2"),
                     labels=c("Acclimation", "Transplant", "Transplant"),
                     values=c("circle","square","square")) +
  geom_point(size=2) +
  geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
  labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
  theme_classic() +
  theme(legend.position="none")

p1<-beta_neutral_nmds_transplant7 %>%
  group_by(type) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
  scale_color_manual(name="Type",
                     breaks=c("Control", "Hot_control", "Treatment"),
                     labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                     values=c("#4477AA","#d57d2c","#76b183")) +
  scale_shape_manual(name="time_point",
                     breaks=c("1_Acclimation", "3_Transplant1", "4_Transplant2"),
                     labels=c("Acclimation", "Transplant", "Transplant"),
                     values=c("circle","square","square")) +
  geom_point(size=2) +
  geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
  labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
  theme_classic() +
  theme(legend.position="none")

p2<-beta_phylo_nmds_transplant7 %>%
  group_by(type) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
  scale_color_manual(name="Type",
                     breaks=c("Control", "Hot_control", "Treatment"),
                     labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                     values=c("#4477AA","#d57d2c","#76b183")) +
  scale_shape_manual(name="time_point",
                     breaks=c("1_Acclimation", "3_Transplant1", "4_Transplant2"),
                     labels=c("Acclimation", "Transplant", "Transplant"),
                     values=c("circle","square","square")) +
  geom_point(size=2) +
  geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
  labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
  theme_classic() +
  theme(legend.position="none")

p3<-beta_func_nmds_transplant7 %>%
  group_by(type) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
  scale_color_manual(name="Type",
                     breaks=c("Control", "Hot_control", "Treatment"),
                     labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                     values=c("#4477AA","#d57d2c","#76b183")) +
  scale_shape_manual(name="time_point",
                     breaks=c("1_Acclimation", "3_Transplant1", "4_Transplant2"),
                     labels=c("Acclimation", "Transplant", "Transplant"),
                     values=c("circle","square","square")) +
  geom_point(size=2) +
  geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
  labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
  theme_classic()+
  theme(legend.position="none")

6.3.5.6 Comparison between Acclimation vs Post-FMT1

post3 <- meta  %>%
  filter(time_point == "1_Acclimation" | time_point == "5_Post-FMT1")

temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

post3.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(post3))]
identical(sort(colnames(post3.counts)),sort(as.character(rownames(post3))))

post3_nmds <- sample_metadata %>%
  filter(time_point == "1_Acclimation" | time_point == "5_Post-FMT1")

6.3.5.7 Number of samples used

[1] 53

beta_div_richness_post3<-hillpair(data=post3.counts, q=0)
beta_div_neutral_post3<-hillpair(data=post3.counts, q=1)
beta_div_phylo_post3<-hillpair(data=post3.counts, q=1, tree=genome_tree)
beta_div_func_post3<-hillpair(data=post3.counts, q=1, dist=dist)

#Arrange of metadata dataframe
post3_arrange<-post3[labels(beta_div_neutral_post3$S),]
post3_arrange$type_time <- interaction(post3_arrange$type, post3_arrange$time_point)

6.3.5.7.1 Richness

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)   
Groups     2 0.099607 0.049803 9.5441    999  0.002 **
Residuals 50 0.260911 0.005218                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
               Control Hot_control Treatment
Control                 0.00300000     0.899
Hot_control 0.00102653                 0.001
Treatment   0.88832670  0.00010131          

	Df	SumOfSqs	R2	F	Pr(>F)
time_point	1	1.2059071	0.0649048	3.913237	0.001
Population	1	1.7615474	0.0948107	5.716321	0.001
time_point:Population	1	0.5122847	0.0275724	1.662393	0.005
Residual	49	15.0998916	0.8127121	NA	NA
Total	52	18.5796308	1.0000000	NA	NA

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	0.3620815	1.052109	0.06169963	0.311	1.000
Control.1_Acclimation vs Hot_control.1_Acclimation	1	1.2800877	4.605444	0.22350616	0.001	0.015	.
Control.1_Acclimation vs Control.5_Post-FMT1	1	0.6845657	1.998114	0.11101796	0.002	0.030	.
Control.1_Acclimation vs Treatment.5_Post-FMT1	1	0.8437461	2.499232	0.14281954	0.002	0.030	.
Control.1_Acclimation vs Hot_control.5_Post-FMT1	1	1.1208022	3.568670	0.18236649	0.001	0.015	.
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	1.3606630	5.087152	0.24124415	0.001	0.015	.
Treatment.1_Acclimation vs Control.5_Post-FMT1	1	0.7216200	2.172734	0.11956009	0.002	0.030	.
Treatment.1_Acclimation vs Treatment.5_Post-FMT1	1	0.9551308	2.926054	0.16322910	0.001	0.015	.
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1	1	1.2263345	4.039487	0.20157637	0.001	0.015	.
Hot_control.1_Acclimation vs Control.5_Post-FMT1	1	1.4319792	5.384836	0.25180628	0.001	0.015	.
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1	1	0.8172413	3.194690	0.17558364	0.001	0.015	.
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.5796135	2.441615	0.13239702	0.001	0.015	.
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1	1	0.5615418	1.729004	0.10335366	0.017	0.255
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.8438429	2.793772	0.14865413	0.002	0.030	.
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.3734921	1.268929	0.07799710	0.111	1.000

6.3.5.7.2 Neutral

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00945 0.0094472 1.1428    999  0.289
Residuals 51 0.42161 0.0082669                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 5_Post-FMT1
1_Acclimation                     0.283
5_Post-FMT1          0.2901            

	Df	SumOfSqs	R2	F	Pr(>F)
time_point	1	1.7277808	0.0986354	6.486824	0.001
Population	1	2.0558578	0.1173647	7.718565	0.001
time_point:Population	1	0.6819354	0.0389303	2.560276	0.004
Residual	49	13.0512643	0.7450696	NA	NA
Total	52	17.5168383	1.0000000	NA	NA

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	0.2316020	0.7712905	0.04598874	0.691	1.000
Control.1_Acclimation vs Hot_control.1_Acclimation	1	1.4015347	5.7562378	0.26457873	0.001	0.015	.
Control.1_Acclimation vs Control.5_Post-FMT1	1	0.8332162	2.9081103	0.15380227	0.001	0.015	.
Control.1_Acclimation vs Treatment.5_Post-FMT1	1	1.1719595	4.0685514	0.21336447	0.002	0.030	.
Control.1_Acclimation vs Hot_control.5_Post-FMT1	1	1.4260875	5.2413171	0.24675104	0.001	0.015	.
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	1.6347704	6.8326887	0.29925029	0.001	0.015	.
Treatment.1_Acclimation vs Control.5_Post-FMT1	1	0.9517634	3.3715700	0.17404733	0.001	0.015	.
Treatment.1_Acclimation vs Treatment.5_Post-FMT1	1	1.3127773	4.6298256	0.23585668	0.001	0.015	.
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1	1	1.6713369	6.2395460	0.28056085	0.001	0.015	.
Hot_control.1_Acclimation vs Control.5_Post-FMT1	1	1.5409781	6.8338056	0.29928456	0.001	0.015	.
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1	1	0.9133614	4.0964534	0.21451383	0.001	0.015	.
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.6954835	3.2951234	0.17077493	0.001	0.015	.
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1	1	0.6051778	2.2508491	0.13047758	0.011	0.165
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	1.0528902	4.1436369	0.20570451	0.001	0.015	.
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.4150076	1.6372683	0.09840968	0.054	0.810

6.3.5.7.3 Phylogenetic

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.05132 0.051320 2.6745    999  0.116
Residuals 51 0.97861 0.019189                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 5_Post-FMT1
1_Acclimation                     0.122
5_Post-FMT1         0.10812            

	Df	SumOfSqs	R2	F	Pr(>F)
time_point	1	0.4329638	0.1243360	7.717455	0.001
Population	1	0.2375991	0.0682323	4.235135	0.003
time_point:Population	1	0.0626513	0.0179918	1.116741	0.270
Residual	49	2.7489923	0.7894398	NA	NA
Total	52	3.4822065	1.0000000	NA	NA

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	0.04186923	0.4391642	0.02671451	0.745	1.000
Control.1_Acclimation vs Hot_control.1_Acclimation	1	0.15609416	2.5546889	0.13768428	0.039	0.585
Control.1_Acclimation vs Control.5_Post-FMT1	1	0.19193367	2.9749922	0.15678490	0.029	0.435
Control.1_Acclimation vs Treatment.5_Post-FMT1	1	0.14627288	1.7907381	0.10665035	0.147	1.000
Control.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.25061348	3.6146185	0.18428187	0.008	0.120
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	0.23108846	4.0521838	0.20208192	0.002	0.030	.
Treatment.1_Acclimation vs Control.5_Post-FMT1	1	0.26358465	4.3608960	0.21417997	0.007	0.105
Treatment.1_Acclimation vs Treatment.5_Post-FMT1	1	0.25319427	3.2738422	0.17915456	0.034	0.510
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.39050120	5.9837393	0.27218933	0.001	0.015	.
Hot_control.1_Acclimation vs Control.5_Post-FMT1	1	0.14203376	5.4200212	0.25303529	0.001	0.015	.
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1	1	0.09666753	2.3682173	0.13635351	0.021	0.315
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.09252600	2.9824958	0.15711821	0.006	0.090
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1	1	0.01842535	0.4144162	0.02688498	0.794	1.000
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.05987967	1.7387847	0.09802164	0.106	1.000
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.03212966	0.6477782	0.04139746	0.691	1.000

6.3.5.7.4 Functional

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00541 0.0054137 0.2021    999  0.659
Residuals 51 1.36615 0.0267873                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 5_Post-FMT1
1_Acclimation                     0.662
5_Post-FMT1         0.65494            

	Df	SumOfSqs	R2	F	Pr(>F)
time_point	1	0.0571984	0.0294799	1.5128525	0.175
Population	1	0.0239890	0.0123639	0.6344904	0.897
time_point:Population	1	0.0064542	0.0033265	0.1707088	0.632
Residual	49	1.8526072	0.9548297	NA	NA
Total	52	1.9402488	1.0000000	NA	NA

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	0.090450145	1.65575459	0.093779882	0.242	1.000
Control.1_Acclimation vs Hot_control.1_Acclimation	1	0.085769225	1.63605364	0.092767559	0.224	1.000
Control.1_Acclimation vs Control.5_Post-FMT1	1	0.030318564	0.53607587	0.032418566	0.529	1.000
Control.1_Acclimation vs Treatment.5_Post-FMT1	1	0.236457683	4.06299320	0.213135113	0.055	0.825
Control.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.135603602	2.22854385	0.122255726	0.160	1.000
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	0.001373886	0.07238159	0.004503476	0.617	1.000
Treatment.1_Acclimation vs Control.5_Post-FMT1	1	0.002173211	0.09402475	0.005842215	0.639	1.000
Treatment.1_Acclimation vs Treatment.5_Post-FMT1	1	0.059119657	2.62461793	0.148917721	0.167	1.000
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.010911935	0.39816986	0.024281360	0.477	1.000
Hot_control.1_Acclimation vs Control.5_Post-FMT1	1	-0.002303582	-0.11016709	-0.006933181	0.753	1.000
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1	1	0.058908140	2.91987617	0.162940644	0.135	1.000
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.005904278	0.23427876	0.014431116	0.543	1.000
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1	1	0.116146557	4.72479132	0.239535681	0.069	1.000
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.050009298	1.70482607	0.096291602	0.229	1.000
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.012358590	0.42381202	0.027477774	0.491	1.000

beta_richness_nmds_post3 <- beta_div_richness_post3$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post3_nmds, by = c("sample" = "Tube_code"))

beta_neutral_nmds_post3 <- beta_div_neutral_post3$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post3_nmds, by = c("sample" = "Tube_code"))

beta_phylo_nmds_post3 <- beta_div_phylo_post3$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post3_nmds, by = join_by(sample == Tube_code))

beta_func_nmds_post3 <- beta_div_func_post3$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post3_nmds, by = join_by(sample == Tube_code))

6.3.5.8 Comparison between Acclimation vs Post-FMT2

post4 <- meta  %>%
  filter(time_point == "1_Acclimation" | time_point == "6_Post-FMT2")

temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

post4.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(post4))]
identical(sort(colnames(post4.counts)),sort(as.character(rownames(post4))))

post4_nmds <- sample_metadata %>%
  filter(time_point == "1_Acclimation" | time_point == "6_Post-FMT2")

6.3.5.9 Number of samples used

[1] 54

beta_div_richness_post4<-hillpair(data=post4.counts, q=0)
beta_div_neutral_post4<-hillpair(data=post4.counts, q=1)
beta_div_phylo_post4<-hillpair(data=post4.counts, q=1, tree=genome_tree)
beta_div_func_post4<-hillpair(data=post4.counts, q=1, dist=dist)

#Arrange of metadata dataframe
post4_arrange<-post4[labels(beta_div_neutral_post4$S),]
post4_arrange$type_time <- interaction(post4_arrange$type, post4_arrange$time_point)

6.3.5.9.1 Richness

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     2 0.06809 0.034047 3.8471    999  0.027 *
Residuals 51 0.45135 0.008850                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              Control Hot_control Treatment
Control                 0.0320000     0.880
Hot_control 0.0349385                 0.005
Treatment   0.8855174   0.0047257          

	Df	SumOfSqs	R2	F	Pr(>F)
time_point	1	0.8124061	0.0462232	2.847438	0.001
Population	1	2.0491994	0.1165926	7.182331	0.001
time_point:Population	1	0.4485668	0.0255219	1.572202	0.003
Residual	50	14.2655595	0.8116623	NA	NA
Total	53	17.5757317	1.0000000	NA	NA

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	0.3620815	1.052109	0.06169963	0.337	1.000
Control.1_Acclimation vs Hot_control.1_Acclimation	1	1.2800877	4.605444	0.22350616	0.001	0.015	.
Control.1_Acclimation vs Treatment.6_Post-FMT2	1	0.8430295	2.845779	0.15100353	0.001	0.015	.
Control.1_Acclimation vs Control.6_Post-FMT2	1	0.5232174	1.683240	0.09518843	0.037	0.555
Control.1_Acclimation vs Hot_control.6_Post-FMT2	1	1.1217138	3.634271	0.18509835	0.001	0.015	.
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	1.3606630	5.087152	0.24124415	0.001	0.015	.
Treatment.1_Acclimation vs Treatment.6_Post-FMT2	1	0.9130048	3.195028	0.16645080	0.001	0.015	.
Treatment.1_Acclimation vs Control.6_Post-FMT2	1	0.5959230	1.984036	0.11032208	0.002	0.030	.
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2	1	1.2747787	4.275366	0.21086503	0.001	0.015	.
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2	1	0.6397330	2.913695	0.15405213	0.001	0.015	.
Hot_control.1_Acclimation vs Control.6_Post-FMT2	1	1.4575447	6.224524	0.28007456	0.001	0.015	.
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2	1	0.3276169	1.412318	0.08111028	0.035	0.525
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2	1	0.6463814	2.560441	0.13795154	0.002	0.030	.
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	0.4796256	1.916520	0.10696943	0.001	0.015	.
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	1.1305044	4.268317	0.21059061	0.001	0.015	.

6.3.5.9.2 Neutral

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.01544 0.0154447 2.0972    999  0.139
Residuals 52 0.38294 0.0073643                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 6_Post-FMT2
1_Acclimation                     0.141
6_Post-FMT2         0.15357            

	Df	SumOfSqs	R2	F	Pr(>F)
time_point	1	1.0151664	0.0602602	3.909554	0.001
Population	1	2.2827471	0.1355037	8.791191	0.001
time_point:Population	1	0.5653146	0.0335570	2.177109	0.001
Residual	50	12.9831505	0.7706790	NA	NA
Total	53	16.8463787	1.0000000	NA	NA

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	0.2316020	0.7712905	0.04598874	0.710	1.000
Control.1_Acclimation vs Hot_control.1_Acclimation	1	1.4015347	5.7562378	0.26457873	0.001	0.015	.
Control.1_Acclimation vs Treatment.6_Post-FMT2	1	1.1746426	4.5564741	0.22165640	0.001	0.015	.
Control.1_Acclimation vs Control.6_Post-FMT2	1	0.5286441	1.9819408	0.11021840	0.003	0.045	.
Control.1_Acclimation vs Hot_control.6_Post-FMT2	1	1.3443224	4.9104417	0.23483204	0.002	0.030	.
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	1.6347704	6.8326887	0.29925029	0.001	0.015	.
Treatment.1_Acclimation vs Treatment.6_Post-FMT2	1	1.3540292	5.3398081	0.25022756	0.001	0.015	.
Treatment.1_Acclimation vs Control.6_Post-FMT2	1	0.6311089	2.4041625	0.13063146	0.002	0.030	.
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2	1	1.6125755	5.9825981	0.27215155	0.001	0.015	.
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2	1	0.6202327	3.1519868	0.16457754	0.001	0.015	.
Hot_control.1_Acclimation vs Control.6_Post-FMT2	1	1.5701179	7.6327037	0.32297209	0.001	0.015	.
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2	1	0.3634438	1.7083388	0.09647087	0.030	0.450
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2	1	1.0227481	4.6483346	0.22511910	0.001	0.015	.
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	0.5010202	2.2065321	0.12119453	0.002	0.030	.
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	1.3619424	5.7710313	0.26507845	0.001	0.015	.

6.3.5.9.3 Phylogenetic

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.06978 0.069777 5.0345    999  0.021 *
Residuals 52 0.72071 0.013860                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 6_Post-FMT2
1_Acclimation                     0.016
6_Post-FMT2        0.029131            

	Df	SumOfSqs	R2	F	Pr(>F)
time_point	1	0.4131659	0.1296152	8.500864	0.001
Population	1	0.2372445	0.0744265	4.881291	0.001
time_point:Population	1	0.1070826	0.0335931	2.203218	0.028
Residual	50	2.4301410	0.7623651	NA	NA
Total	53	3.1876340	1.0000000	NA	NA

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	0.04186923	0.4391642	0.02671451	0.723	1.000
Control.1_Acclimation vs Hot_control.1_Acclimation	1	0.15609416	2.5546889	0.13768428	0.035	0.525
Control.1_Acclimation vs Treatment.6_Post-FMT2	1	0.26322331	4.3060281	0.21205664	0.004	0.060
Control.1_Acclimation vs Control.6_Post-FMT2	1	0.16047895	2.5405742	0.13702781	0.045	0.675
Control.1_Acclimation vs Hot_control.6_Post-FMT2	1	0.25529510	4.0109138	0.20043631	0.002	0.030	.
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	0.23108846	4.0521838	0.20208192	0.003	0.045	.
Treatment.1_Acclimation vs Treatment.6_Post-FMT2	1	0.36496892	6.3966666	0.28560797	0.002	0.030	.
Treatment.1_Acclimation vs Control.6_Post-FMT2	1	0.22628210	3.8292220	0.19311005	0.023	0.345
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2	1	0.34830814	5.8463335	0.26761166	0.003	0.045	.
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2	1	0.10002871	4.3836237	0.21505615	0.001	0.015	.
Hot_control.1_Acclimation vs Control.6_Post-FMT2	1	0.12577510	5.0601287	0.24027055	0.001	0.015	.
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2	1	0.06334378	2.4997737	0.13512455	0.031	0.465
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2	1	0.05927454	2.3820253	0.12958449	0.018	0.270
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	0.06906280	2.7224602	0.14541146	0.005	0.075
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	0.11081709	4.0436561	0.20174244	0.001	0.015	.

6.3.5.9.4 Functional

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00538 0.0053768 0.1915    999  0.681
Residuals 52 1.46001 0.0280772                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 6_Post-FMT2
1_Acclimation                     0.675
6_Post-FMT2         0.66348            

	Df	SumOfSqs	R2	F	Pr(>F)
time_point	1	0.0451272	0.0242533	1.2655291	0.234
Population	1	0.0004181	0.0002247	0.0117261	0.288
time_point:Population	1	0.0321822	0.0172961	0.9025046	0.294
Residual	50	1.7829366	0.9582260	NA	NA
Total	53	1.8606640	1.0000000	NA	NA

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	9.045014e-02	1.65575459	0.0937798825	0.236	1
Control.1_Acclimation vs Hot_control.1_Acclimation	1	8.576922e-02	1.63605364	0.0927675587	0.243	1
Control.1_Acclimation vs Treatment.6_Post-FMT2	1	1.171666e-01	2.16526354	0.1191980254	0.172	1
Control.1_Acclimation vs Control.6_Post-FMT2	1	1.134013e-01	2.16495643	0.1191831338	0.178	1
Control.1_Acclimation vs Hot_control.6_Post-FMT2	1	6.533744e-02	0.94814395	0.0559438221	0.272	1
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	1.373886e-03	0.07238159	0.0045034760	0.649	1
Treatment.1_Acclimation vs Treatment.6_Post-FMT2	1	-4.929805e-04	-0.02385162	-0.0014929516	0.697	1
Treatment.1_Acclimation vs Control.6_Post-FMT2	1	2.443617e-03	0.12903840	0.0080003779	0.538	1
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2	1	6.968380e-03	0.19647180	0.0121305310	0.562	1
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2	1	3.866607e-04	0.02093980	0.0013070267	0.652	1
Hot_control.1_Acclimation vs Control.6_Post-FMT2	1	-5.633463e-05	-0.00336651	-0.0002104511	0.691	1
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2	1	2.011448e-03	0.06046867	0.0037650628	0.727	1
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2	1	-8.527330e-03	-0.46290555	-0.0297935723	0.836	1
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	-1.648721e-03	-0.04717131	-0.0029569243	0.905	1
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	4.367477e-03	0.13147026	0.0081499244	0.676	1

beta_richness_nmds_post4 <- beta_div_richness_post4$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post4_nmds, by = c("sample" = "Tube_code"))

beta_neutral_nmds_post4 <- beta_div_neutral_post4$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post4_nmds, by = c("sample" = "Tube_code"))

beta_phylo_nmds_post4 <- beta_div_phylo_post4$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post4_nmds, by = join_by(sample == Tube_code))

beta_func_nmds_post4 <- beta_div_func_post4$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post4_nmds, by = join_by(sample == Tube_code))

6.3.6 6. Are there differences between the control and the treatment group?

6.3.6.1 After 1 week –> Post-FMT1

post1 <- meta %>%
  filter(time_point == "5_Post-FMT1")

post1.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(post1))]
identical(sort(colnames(post1.counts)),sort(as.character(rownames(post1))))

post1_nmds <- sample_metadata %>%
  filter(time_point == "5_Post-FMT1")

6.3.6.2 Number of samples used

[1] 26

beta_div_richness_post1<-hillpair(data=post1.counts, q=0)
beta_div_neutral_post1<-hillpair(data=post1.counts, q=1)
beta_div_phylo_post1<-hillpair(data=post1.counts, q=1, tree=genome_tree)
beta_div_func_post1<-hillpair(data=post1.counts, q=1, dist=dist)

#Arrange of metadata dataframe
post1_arrange<-post1[labels(beta_div_neutral_post1$S),]

6.3.6.2.1 Richness

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.017675 0.0088373 2.3825    999  0.103
Residuals 23 0.085312 0.0037092                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              Control Hot_control Treatment
Control                 0.0030000     0.665
Hot_control 0.0068795                 0.216
Treatment   0.6248469   0.2084296          

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	0.6340254	0.0768024	2.065607	0.005
type	1	0.5615418	0.0680222	1.829461	0.008
Residual	23	7.0597099	0.8551754	NA	NA
Total	25	8.2552771	1.0000000	NA	NA

                     pairs Df SumsOfSqs  F.Model        R2 p.value p.adjusted sig
1     Control vs Treatment  1 0.5615418 1.729004 0.1033537   0.014      0.042   .
2   Control vs Hot_control  1 0.8438429 2.793772 0.1486541   0.001      0.003   *
3 Treatment vs Hot_control  1 0.3734921 1.268929 0.0779971   0.112      0.336    

6.3.6.2.2 Neutral

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.011001 0.0055005 0.6303    999  0.569
Residuals 23 0.200714 0.0087267                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.22800     0.960
Hot_control 0.21166                 0.452
Treatment   0.95468     0.43604          

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	0.7907904	0.1076445	3.056657	0.001
type	1	0.6051778	0.0823784	2.339205	0.006
Residual	23	5.9503501	0.8099772	NA	NA
Total	25	7.3463184	1.0000000	NA	NA

                     pairs Df SumsOfSqs  F.Model         R2 p.value p.adjusted sig
1     Control vs Treatment  1 0.6051778 2.250849 0.13047758   0.017      0.051    
2   Control vs Hot_control  1 1.0528902 4.143637 0.20570451   0.003      0.009   *
3 Treatment vs Hot_control  1 0.4150076 1.637268 0.09840968   0.056      0.168    

6.3.6.2.3 Phylogenetic

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.00440 0.0021994 0.1369    999  0.912
Residuals 23 0.36941 0.0160614                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.94100     0.686
Hot_control 0.91505                 0.756
Treatment   0.63312     0.73046          

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	0.0560850	0.0531376	1.3149967	0.238
type	1	0.0184254	0.0174571	0.4320099	0.810
Residual	23	0.9809570	0.9294053	NA	NA
Total	25	1.0554673	1.0000000	NA	NA

                     pairs Df  SumsOfSqs   F.Model         R2 p.value p.adjusted sig
1     Control vs Treatment  1 0.01842535 0.4144162 0.02688498   0.775      1.000    
2   Control vs Hot_control  1 0.05987967 1.7387847 0.09802164   0.128      0.384    
3 Treatment vs Hot_control  1 0.03212966 0.6477782 0.04139746   0.690      1.000    

6.3.6.2.4 Functional

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq     F N.Perm Pr(>F)
Groups     2 0.00400 0.0020014 0.145    999   0.87
Residuals 23 0.31753 0.0138057                    

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.60200     0.746
Hot_control 0.59817                 0.840
Treatment   0.75141     0.83718          

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	0.0024979	0.0033024	0.0900845	0.608
type	1	0.1161466	0.1535542	4.1887855	0.087
Residual	23	0.6377435	0.8431434	NA	NA
Total	25	0.7563879	1.0000000	NA	NA

                     pairs Df  SumsOfSqs  F.Model         R2 p.value p.adjusted sig
1     Control vs Treatment  1 0.11614656 4.724791 0.23953568   0.071      0.213    
2   Control vs Hot_control  1 0.05000930 1.704826 0.09629160   0.217      0.651    
3 Treatment vs Hot_control  1 0.01235859 0.423812 0.02747777   0.453      1.000    

beta_richness_nmds_post1 <- beta_div_richness_post1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post1_nmds, by = join_by(sample == Tube_code))

beta_neutral_nmds_post1 <- beta_div_neutral_post1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post1_nmds, by = join_by(sample == Tube_code))

beta_phylogenetic_nmds_post1 <- beta_div_phylo_post1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post1_nmds, by = join_by(sample == Tube_code))

beta_functional_nmds_post1 <- beta_div_func_post1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post1_nmds, by = join_by(sample == Tube_code))

p0<-beta_richness_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p1<-beta_neutral_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")
  
p2<-beta_phylogenetic_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p3<-beta_functional_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
                theme_classic()+
                theme(legend.position="none")

ggarrange(p0, p1, p2, p3, ncol=2, nrow=2, common.legend = TRUE, legend="right")

6.3.6.3 After 2 weeks –>Post-FMT2

post2 <- meta %>%
  filter(time_point == "6_Post-FMT2")

post2.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(post2))]
identical(sort(colnames(post2.counts)),sort(as.character(rownames(post2))))

post2_nmds <- sample_metadata %>%
  filter(time_point == "6_Post-FMT2")

6.3.6.4 Number of samples used

[1] 27

beta_div_richness_post2<-hillpair(data=post2.counts, q=0)
beta_div_neutral_post2<-hillpair(data=post2.counts, q=1)
beta_div_phylo_post2<-hillpair(data=post2.counts, q=1, tree=genome_tree)
beta_div_func_post2<-hillpair(data=post2.counts, q=1, dist=dist)

#Arrange of metadata dataframe
post2_arrange<-post2[labels(beta_div_neutral_post2$S),]

6.3.6.4.1 Richness

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.002011 0.0010056 0.1982    999  0.823
Residuals 24 0.121775 0.0050740                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.70400     0.826
Hot_control 0.67789                 0.609
Treatment   0.79246     0.59820          

	Df	SumOfSqs	R2	F	Pr(>F)
type	2	1.504341	0.1967776	2.939822	0.001
Residual	24	6.140538	0.8032224	NA	NA
Total	26	7.644879	1.0000000	NA	NA

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Treatment vs Control	1	0.6463814	2.560441	0.1379515	0.002	0.006	*
Treatment vs Hot_control	1	0.4796256	1.916520	0.1069694	0.001	0.003	*
Control vs Hot_control	1	1.1305044	4.268317	0.2105906	0.001	0.003	*

6.3.6.4.2 Neutral

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.008262 0.0041311 0.8024    999  0.477
Residuals 24 0.123559 0.0051483                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.45800     0.643
Hot_control 0.44675                 0.273
Treatment   0.65989     0.25095          

	Df	SumOfSqs	R2	F	Pr(>F)
type	2	1.923807	0.2603795	4.224537	0.001
Residual	24	5.464666	0.7396205	NA	NA
Total	26	7.388473	1.0000000	NA	NA

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Treatment vs Control	1	1.0227481	4.648335	0.2251191	0.001	0.003	*
Treatment vs Hot_control	1	0.5010202	2.206532	0.1211945	0.002	0.006	*
Control vs Hot_control	1	1.3619424	5.771031	0.2650785	0.001	0.003	*

6.3.6.4.3 Phylogenetic

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.000407 0.0002034 0.0487    999  0.963
Residuals 24 0.100305 0.0041794                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.93200     0.845
Hot_control 0.93765                 0.756
Treatment   0.83933     0.76015          

	Df	SumOfSqs	R2	F	Pr(>F)
type	2	0.1594363	0.2042241	3.079623	0.001
Residual	24	0.6212564	0.7957759	NA	NA
Total	26	0.7806927	1.0000000	NA	NA

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Treatment vs Control	1	0.05927454	2.382025	0.1295845	0.022	0.066
Treatment vs Hot_control	1	0.06906280	2.722460	0.1454115	0.004	0.012	.
Control vs Hot_control	1	0.11081709	4.043656	0.2017424	0.002	0.006	*

6.3.6.4.4 Functional

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.01126 0.0056302 0.2861    999  0.814
Residuals 24 0.47233 0.0196806                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.53800     0.671
Hot_control 0.48255                 0.812
Treatment   0.60116     0.75643          

	Df	SumOfSqs	R2	F	Pr(>F)
type	2	-0.0038724	-0.0056213	-0.0670788	0.921
Residual	24	0.6927468	1.0056213	NA	NA
Total	26	0.6888744	1.0000000	NA	NA

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Treatment vs Control	1	-0.008527330	-0.46290555	-0.029793572	0.873	1
Treatment vs Hot_control	1	-0.001648721	-0.04717131	-0.002956924	0.899	1
Control vs Hot_control	1	0.004367477	0.13147026	0.008149924	0.669	1

beta_richness_nmds_post2 <- beta_div_richness_post2$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post2_nmds, by = join_by(sample == Tube_code))

beta_neutral_nmds_post2 <- beta_div_neutral_post2$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post2_nmds, by = join_by(sample == Tube_code))

beta_phylogenetic_nmds_post2 <- beta_div_phylo_post2$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post2_nmds, by = join_by(sample == Tube_code))

beta_functional_nmds_post2 <- beta_div_func_post2$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post2_nmds, by = join_by(sample == Tube_code))

p0<-beta_richness_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p1<-beta_neutral_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")
  
p2<-beta_phylogenetic_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p3<-beta_functional_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
                theme_classic()+
                theme(legend.position="none")

ggarrange(p0, p1, p2, p3, ncol=2, nrow=2, common.legend = TRUE, legend="right")